读书笔记《Matlab在语音信号分析与合成中的应用》第5章带噪语音和预处理

5.1 纯语音和带噪语音

噪声的分类
- 加性噪声
  - $x(n)=s(n)+d(n)$
- 非加性噪声
  - 卷机性噪声
  - 乘积性噪声

5.2 信噪比

公式

$SNR=10log_{10}\frac{\sum_{n=0}^{N-1}s^2(n)}{\sum_{n=0}^{N-1}d^2(n)}$

信号的能量比上噪声的能量，再求分贝

5.3 Matlab产生带噪语音

构造函数Gnoisegen

function [y,noise] = Gnoisegen(x,snr)
#1. 利用随机函数构造noise
noise_pre = randn(size(x));
#2. 为了便于测试，需要设定一个参数SNR，计算达到SNR要求的noise的能量
Nx = length(x);
signal_power = 1/Nx*sum(x.*x);
noise_pre_power = 1/Nx*sum(noise_pre.*noise_pre);
noise_power = signal_power/(10^(snr/10));
noise_ampli_cof = sqrt(noise_power/noise_pre_power);
#3. 将构造好的noise加到原信号上
y = x + noise_ampli_cof * noise;

计算SNR，需要传入原始信号与带噪信号

function snr = SNR_singlech(signal,s_plus_n)
signal_power = sum((signal-mean(signal).^2));
noise_power = sum((s_plus_n-signal).^2);
# 分贝的公式
snr = 10 * log10(signal_power/noise_power);

将手中有的任意的noise，叠加到原信号中：

采样频率不同
长度不同

function [signal,noise]=add_noisedata(s,data,fs,fs1,snr)
s=s(:);                          % 把信号转换成列数据
s=s-mean(s);                     % 消除直流分量
sL=length(s);                    % 求出信号的长度

if fs~=fs1                       % 若纯语音信号的采样频率与噪声的采样频率不相等
    x=resample(data,fs,fs1);     % 对噪声重采样，使噪声采样频率与纯语音信号的采样频率相同
else
    x=data;
end

x=x(:);                          % 把噪声数据转换成列数据
x=x-mean(x);                     % 消除直流分量
xL=length(x);                    % 求噪声数据长度
if xL>=sL                        % 如果噪声数据长度与信号数据长度不等，把噪声数据截断或补足
    x=x(1:sL);
else
    disp('Warning: 噪声数据短于信号数据，以补0来补足！')
    x=[x; zeros(sL-xL,1)];
end

Sr=snr;
Es=sum(x.*x);                    % 求出信号的能量
Ev=sum(s.*s);                    % 求出噪声的能量
a=sqrt(Ev/Es/(10^(Sr/10)));      % 计算出噪声的比例因子
noise=a*x;                       % 调整噪声的幅值
signal=s+noise;                  % 构成带噪语音

5.4 语音信号的预处理–多项式拟合消除趋势项

消除趋势项和直流分量
- 最小二乘法拟合趋势项【这里全部用的是解析法去让二次项最小】
  - detrend函数 - 消除线性趋势项
    
    $y_k = x_k - \hat{x_k}=x_k - (a_0+a_1*k)$
  - polydetrend函数
    [y,xtrend] = polydetrend(x,fs,m)
    # xtrend 是叠加在序列上的趋势项序列。
    # m 可以控制拟合的次数
    $\hat{x}_{k}=a_{0}+a_{1} k+a_{2} k^{2}+\cdots+a_{m} k^{m}=\sum_{j=0}^{m} a_{j} k^{j} \quad k=1,2, \cdots, n$

将拟合后的xtrend减去就得到 $y$

$y_k = x_k - \hat{x_k}$

其本质原理就是两步：

调用polyfit最小二乘法拟合出 $\hat{x_t}$ 的系数a
调用polyval(a,x)构造出 $\hat{x_t}$ 序列，即 $xtrend$ 趋势项序列

原信号减去趋势项序列得到y

function [y,xtrend]=polydetrend(x, fs, m)
x=x(:);                 % 把语音信号x转换为列数据
N=length(x);            % 求出x的长度
t= (0: N-1)'/fs;        % 按x的长度和采样频率设置时间序列
a=polyfit(t, x, m);     % 用最小二乘法拟合语音信号x的多项式系数a
xtrend=polyval(a, t);   % 用系数a和时间序列t构成趋势项
y=x-xtrend;             % 从语音信号x中清除趋势项
end

5.5 语音信号的预处理–滤波器的设计

基音频率范围大部分在60~450Hz范围之间，在采集到语音信号后为了提取基音只需要语音的低频成分即可，故需要低频滤波器，同时，采样过程中非常容易把公频50Hz的交流声给混入到语音信号中，故还需要高通滤波器除50Hz的干扰

IIR滤波器
- 由模拟滤波器的原型变换而来
  - 巴特沃斯滤波器
  - 切比雪夫I型
  - 切比雪夫II型
  - 椭圆型滤波器
- 具体使用方法
  1. 本质上就是确认滤波器的系数 a和b
    
    $H(z)=\frac{B(z)}{A(z)}=\frac{b(1)+b(2) z^{-1}+\cdots+b(n+1) z^{-n}}{1+a(2) z^{-1}+\cdots+a(n+1) z^{-n}}$
  2. 上面的系数可由滤波器相关函数得到
    [b,a]=cheby2(n,Rsm,Wn) # 切比雪夫II型滤波器
    y = filter(b,a,signal)
FIR 滤波器
- 核心公式
  
  $h(n)=h_d(n)W_d(n)$
  
  h(n)是滤波器的单位脉冲响应
  
  $h_d(n)$ 理想频率响应下求出的单位脉冲响应
  
  $W_d(n)$ 是窗函数

窗函数

Matlab中窗函数的类型

Wd = boxcar(N) %矩形窗
Wd = triang(N) %三角窗
Wd = hannig(N) %海宁窗
Wd = hamming(N) %汉明窗
Wd = blackman(N) %布莱克曼窗
Wd = kaiser(N,beta) %凯泽窗

窗函数的选用依据

窗函数	旁瓣/dB	近似过渡带宽	精确过渡带宽	阻带最小衰减/dB
矩形窗	-13	4 $\pi/N$	1.8 $\pi/N$	21
三角形窗	-25	8 $\pi/N$	6.1 $\pi/N$	25
海宁窗	-31	8 $\pi/N$	6.2 $\pi/N$	44
汉明窗	-41	8 $\pi/N$	6.6 $\pi/N$	53
布莱克曼窗	-57	12 $\pi/N$	11 $\pi/N$	74
泽凯窗 $(\beta=7.865)$	-57		10 $\pi/N$	80